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Space and time discretisations of parabolic differential equations with dynamic boundary conditions are 
studied in a weak formulation that fits into the standard abstract formulation of parabolic problems, 
just that the usual L,2(£2) inner product is replaced by an L2(£2) ®L2(d£2) inner product. The class of 
parabolic equations considered includes linear problems with time- and space-dependent coefficients and 
semi-linear problems such as reaction-diffusion on a surface coupled to diffusion in the bulk. The spatial 
discretisation by finite elements is studied in the proposed framework, with particular attention to the 
error analysis of the Ritz map for the elliptic bilinear form in relation to the inner product, both of which 
contain boundary integrals. The error analysis is done for both polygonal and smooth domains. We 
further consider mass lumping, which enables us to use exponential integrators and bulk-surface splitting 
for time integration. 
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1. Introduction 


We are interested in the numerical solution of parabolic initial-boundary value problems with dynamic 
boundary conditions. Prototypes for this class of problems are the heat equation with Wentzell boundary 
conditions 


{ d t u = Au in £2 

I p d,u = — Ku — d v u on T, 


( 1 . 1 ) 


set on a bounded, piecewise smooth domain £2 C W 1 (the bulk) with boundary T = d£2 (referred to as 
the surface), with a positive coefficient /i and real K and with d v u denoting the normal derivative of u 
on r ; and diffusion on the surface coupled to diffusion in the bulk , 


d t u = Au in £2 

p,d t u = /3 Apu — d v ii on r, 


( 1 . 2 ) 
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with positive coefficients /J. and /3 and with the Laplace-Beltrami operator Ap. We will study numerical 
methods for such equations as well as for inhomogeneous, non-autonomous and semi-linear variants. 

It turns out that such problems admit a weak formulation that fits into the standard abstract frame¬ 
work of parabolic problems. The difference to the heat equation with homogeneous Neumann boundary 
conditions is mainly that the role of the Hilbert spaces Ljifl) and II 1 (£2 j in the Neumann problem is 
taken by other Hilbert spaces, for the above problems by L2(£2) CD LiiT j and an appropriate subspace 
ofH\£2), respectively. 

The papers by Cherfils et al. ( 2010|) and Cherfils & Petcul (|2014 ) on the finite element discretisa¬ 
tion of the Cahn-Hilliard equation on a slab with various dynamic boundary conditions are the only 
publications on the numerical analysis of parabolic problems with dynamic boundary conditions that 
we know of. On the other hand, there has been much r ecent work o n analy t ical and model l ing as ¬ 
pects of such prob lems; s ee, e.g., Cavaterra et al. d2.0 1 Oh : ICoclite et al} d20(>4: IColli & Fukaol (|2014 |) : 
Engel & Fraenellil (l2005l): Favini et al. ( 2002h: Galf i 2008 ): Gal & Grassellil ( 2008 ): Goldstein ( 2006 ) 


Goldstein et al. ( 2011 ): Kenzler et al. ( 200 ll) : Fiero ( 2013 ): Racke & Zheng 12003 ): Vazquez & Vitillaro 

d201lh . We are not aware that the abstract framework of the present paper is common in the literature, 
but related variational settings do appear in some of the references. 

In Section[2]we discuss the abstract variational framework and how problems such as (11.11 ) and (11.21) 
fit in, as well as variants of these problems with space- and time-dependent coefficient functions and 
semi-linear problems such as reaction-diffusion on the boundary coupled to diffusion in the bulk, the 
Allen-Cahn equation with dynamic boundary conditions, and the Cahn-Hilliard equation with various 
types of dynamic boundary conditions. 

In Sections 3 and 4 we study the finite element semi-discretisation in space, which for problems 
(11.11) and (11.21) leads to a large system of ordinary differential equations for the nodal vector u(f), 


Mu(f) + Au(f) = 0, 


where the positive definite mass matrix M and the positive definite or semi-definite stiffness matrix A 
differ from those of the heat equation with Neumann boundary conditions only in entries that correspond 
to pairs of adjacent boundary nodes. Much of the standard numerical analysis of the heat equation with 
Dirichlet or Neumann boundary conditions carries over in a direct way, since the problems with dynamic 
boundary conditions share the same abstract variational framework. There are, however, a few issues 
where care is needed in the extension of the theory: 

• The Ritz projection is based on a different elliptic form that contains boundary integrals. It must 
be put in relation with a different inner product that also contains boundary integrals. 

• Mass lumping is done for an inner product with boundary terms. It behaves differently for ( 11.11) 
than for (11.21) . 

• In the case of non-polygonal domains, the boundary approximation may play a more important 
role than for pure bulk problems. Its effect on the approximation needs to be studied thoroughly. 

We first study a class of linear constant-coefficient problems that includes (11.11 ) and (11.21) . then turn 
to problems with space- and time-dependent coefficient functions and finally to semi-linear problems 
such as the Allen-Cahn equation with dynamic boundary conditions. We discuss finite element space 
discretisation for the case of polygonal domains in Section 3 and then extend the results to smooth 
domains in Section 4. 
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In Section 5 we turn to time discretisation. Known stability and approximation results for standard 
implicit integrators such as backward difference formulae and Radau IIA implicit Runge-Kutta methods 
extend from Dirichlet or Neumann boundary conditions to dynamic boundary conditions without much 
ado, again thanks to the common abstract framework. This framework also makes it obvious how to ap¬ 
ply exponential integrators to the problem class studied here, which is not immediately evident from the 
strong formulation (II.Il l or (11.21 ). We further consider two classes of bulk-surface splitting integrators, 
force splitting and component splitting, where in both cases differential equations corresponding to the 
interior domain and to the boundary are solved separately in an alternating way. 


2. Variational formulation of parabolic problems with dynamic boundary conditions 

2.1 Linear problems with time-independent operators 

2.1.1 Abstract setting. We recall the usual weak formulation of abstract linear parabolic problems: 
Given are two Hilbert spaces V, with norm || • ||, and H, with norm | • | corresponding to the inner product 
(•,•), such that V is densely and continuously embedded in H. On V we consider a continuous bilinear 
form a(-, •) that satisfies a Garding inequality: there exist a > 0 and real c such that 

a(v,v) ^ a||v|| 2 — c|v| 2 Vv£ V. 


On a time interval 0 ^ t ^ T, for given initial data wo £ // and an inhomogeneity / £ LjfO.T\H), the 
abstract parabolic initial value problem then reads: Find u £ C( [(),/'],//) 01^(0, T;V) such that (with 
ii = du/dt) 

{u(t),v) + a(u(t),v) = (f(t),v) VveV (0<f<r) 

u( 0) = Mo • 


We note that (for / = 0) this can be viewed as the //-gradient flow, (w, v) = —E'(u)v for all v £ V, for 
the quadratic energy functional 


E(v) = ja(v,v), vGV. 

The well-posedness of this abstract problem is well known; see, e.g., Dautrav & Lionsl ( 1992 ): Katol 
dl995h . An a priori estimate of the solution is obtained by the familiar energy technique: test with 
v = u(t), note ( u{t),u{t )) = \^W(t)\ 2 and integrate fromO to t to obtain 


j\u(t)\ 2 - ±\u 0 \ 2 + 2 j E(u(s))ds = J (f(s),u(s))ds. 


Using here the bound |(/(s),m(s))| ^ ||/(s)||* ||m(s)|| with the dual norm ||<p||* = sup|| v || =1 |(<p,v)| and 
estimating this further as |(/(s),m(s))| ^ ^||/(s)||* + ■§ll M (' s )l | 2 an d using the Garding inequality, one 
arrives at 

\u(t)\ 2 + a [ t \\u(s)\\ 2 ds^\u 0 \ 2 + ^ [ t \\f(s)\\Us + 2c f \u(s)\ 2 ds, 

Jo a, Jo Jo 

and Gronwall’s inequality then yields 


l“( ? )l 2 + <X L ll M ( s )ll 2ds ^ e2C, (l“°l 2 + lcJ 0 (2 - 2 ^ 

The standard example is V = // ( | (12) and // = L? (G) with the Dirichlet energy functional E(v) = 
ill V l ’ll/ / ,(_o)' i n which case ( 12.11 ) yields the weak formulation of the heat equation on the domain £2 with 
homogeneous Dirichlet boundary conditions. 















4 


B. KOVACS AND CH. LUBICH 


2.1.2 Weak formulation of the heat equation with dynamic boundary conditions. As it turns out, 
the weak formulation of the heat equation with dynamic boundary conditions (II. Il l and (11.2b fits into 
the same abstract framework with the particular choice of Hilbert space V = {v G // 1 (22 ) : sjp yv G 
H l (T)}, where yv denotes the trace of v on the boundary F = r)Q. and with the bilinear form on V 

a(u,v) = f Vw-Vvcbc+tc [ (yu)(yv)da + P I Vr(yu) • Vr(yv)d<7 (kg K,/3>0), (2.3) 

Jn Jr Jr 

where V r (yy) = (I — w r )y(Vv) is the tangential gradient on F (v denotes the unit normal on F), which 
is known to depend only on the trace yv on F. For brevity we will write V r v instead of V r (yv) in the 
following. We have j3 = 0 for (11.11) and K 0 for (11.21 ). 

The Hilbert space H is the completion of V with respect to the norm induced by the inner product 

(n,v) = / uvdLx + p / (■ yu)(yv)da (ju. > 0), (2.4) 

Jn Jr 

so that H is (isomorphic to) Lo(Q)®L 2 (r). 

The corresponding energy functional is then 

E ( v ) = 2 ll Vv ’Hi(i 2 ) + 2 K 'll7'’llL(r) + 3^ll V ^ l, llL(r)> and hf = IMIl^aj + MINI l 2 (ry 

We thus obtain the energy estimate, here stated for the weak solution u(t) = u(-,t ) of (12.1b with / = 0 
for simplicity, 

IWOllL^) + Mll7M(0llL(r) + J Q (ll v «( s )ll l 2 (n) + K \\Ms)\\l 2 (r)+P ll v rM(s)llL ( r)) d? 

^ ll M (0)llz^fra) d Ml!T M (0)|| L 2 (r)- 

The relationship with the strong formulation of the heat equation with dynamic boundary conditions 
(II.lb and (11.2b is given by the following result. 

Lemma 2.1 Every classical solution u G C 2 (£2 x [0,7’]) of the heat equation with dynamic boundary 
conditions 

{ d t u = Au in Q. 

fj.d t u = —Ku + P Apu — d v u on F 

with initial data no is a solution of the weak formulation ( 12.1b with the bilinear forms a(-,-) and (■,•) 
of (12.3b —( 12.4b . Conversely, if the solution n of (12.1b with a(-. ■) and (•,■) of ( 12.3b — (12.4b is sufficiently 
regular, then n is a solution of the strong formulation ( 12.5b . 

The proof is almost identical to the proof of the corresponding result for the heat equation with 
Neumann boundary conditions. It is again based on Green’s formula in the domain and on the boundary, 
and on the fundamental lemma of variational calculus. The proof is therefore omitted. 

2.2 Linear problems with time-varying operators 

2.2.1 Abstract setting. Consider again Hilbert spaces V and H such that V is continuously and 
densely embedded in H. On V we consider uniformly equivalent time-dependent inner product norms 
|| • || f , and on H a time-dependent family of inner products m(t\ ■, ■): H x H —> R that induce uniformly 
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equivalent norms \w\^ = m(t',w,w) for w £H. We assume a bounded partial derivative of m with respect 
to t: 

I dm I 

\ — (f,u,v)\^M' 0 \u\ t \v\ t Vu,veV (0 sSfsSF). 


dt 

On V we consider a time-dependent family of bilinear forms a(t;-,-) : V x V - 
satisfy a uniform Garding inequality 

a(f,v,v) > a||v||^ — c|v|^ Vv€V (0<f^F) 

and are uniformly bounded, 

|a(/;n,v)| < Mi||n|| r ||v|| f Vn,vSV (0 

We further assume 

^(t;u,v) <Afi ||n||, ||v|| f Vn,v€V (0 


(2.6) 

(0 ^ t ^ T) that 


(2.7) 


The non-autonomous version of the abstract parabolic initial value problem is then to find u £ C( [0, T] ,H) 0 
Li (0, T; V) such that 


m(f,u(t),v)+a(t;u(t),v) = m(t;f(t),v) VveV (0 <t^T) 
u( 0) = u o. 


(2.8) 


Here again we obtain an a priori bound using the energy technique: by the same arguments as in the 
time-invariant situation and using the above bounds we obtain 


Ht)\? + aJ ||M(s)||?ds<e 2(c+M o)^| Mo |o + ^ ll/WllL d? 

with the time-dependent dual norm ||<p||* tf = sup|| v || f=1 \m(t; (p, v)|. 


(2.9) 


2.2.2 Heat equation with non-autonomous dynamic boundary conditions. We consider the weak for¬ 
mulation of the problem 


d t u = Au 

p(x,t)d t u(x,t) = — K(x,t)u(x,t) + Vj- • j3 (x,t)Vru{x,t) ~ d v u(x,t) 


in Q 

on r = dQ. 


( 2 . 10 ) 


with real-valued coefficient functions p, K, ft onT x [0, T] such that p,K,fi : [0, T] —>Loo(M) are continu¬ 
ously differentiable, p has a strictly positive lower bound, and p either has a strictly positive lower bound 
or vanishes identically. This fits into the above framework for V = {v £ H 1 (£2) : \//3 Vj-(yv) £Lo(r)} 
with 

a(t;u,v) = [ Vn-Vvdr + f (yM)(yv)d<7 + [ Vpu- Vj-vda (2.11) 

Jn Jr Jr 

and H = L 2 (Q) (BL 2 (T) with 

m(t;u,v)= / mv dr T / p(-,t) (yu)(yv)da. (2.12) 

Jn Jr 

The framework applies equally when p (•, t ) is piecewise continuous on F and has a positive lower bound 
on a subset r + C F and is zero on the complementary part To = F\F|_. In this case, H = La{Q)®L 2 {r+). 
Similarly, ft may be allowed to have a positive lower bound on a time-independent subset of F and to 
vanish on the complement. 
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2.3 Some nonlinear examples 

We present three examples of nonlinear parabolic equations with dynamic boundary conditions. They 
can all be cast in the following abstract form of a semi-linear parabolic problem on suitable spaces V 
and H: Find u £ C([0,T],H) (1L2(0, T\V) such that (with u = du/dt) 

(u(t),v) + a(u(t),v) = {f(u(t)),v) VvGV (0 <t^T) 

U{ 0) = WQ , 


where / : V —► H is a sufficiently regular nonlinearity. In our examples the bilinear form a(-, ■) on V and 
the inner product (-, •) on H contain boundary terms. The bilinear form at-. •) is of the form (12.3b with 
ft > 0. In the first two examples the inner product (•, •) is of the form (12.4b on H = L 2 (Q) ® LjfT), and 
in the third example on H = ©Li(T) or H = 


2.3.1 Reaction-diffusion on a surface coupled to diffusion in the bulk. We consider a reaction-diffusion 
equation on the boundary 

jU<9,y/ = j3 4 r v/ + /(y/) onT, 

with p ,/3>0 and a nonlinear pointwise reaction term f(tj))(x) = f{(j>{x)) (for a smooth and bounded 
function / : M —> R, say), and we further consider diffusion in the bulk, 

d t u = Au inf2. 

We couple these equations subject to the constraint y/ = yu. We obtain the following weak formulation: 
find (m, y/j : [0,T] —> H 1 (£2) x H 1 (F) subject to y/ = yu such that for all (v,0) €H 1 (£2) xH l (r) with 
0 = yv. 


(“> v )z. 2 (t2) +M(V / 5 0)L 2 (r) — -(VM,Vv) L2 (i 2 ) - p(V r yr, V r 0) i2 ( r ) + {f{y>)A)L 2 (r)- 

Equivalently, with the bilinear forms 

(m,v) = {u,v) L2in) +p(yu,y\’) L2(r) 4 

a( u ,v) = (Vn,Vv) L2 (i 2 )+j3(V r M,V r v) L2 (r) 

on the Hilbert spaces // = L^JQ) ©Lt(F) and V = {v £ // 1 (X2) : yv £ // 1 (X)}, respectively, we have 


(u,v) + a(u,v) = (f(yu),yv) L2(n Vv£ V. 


This is to be solved for u £ C(\(),T],Ii) (TX(0, T;V) for given initial data uq £ H. The corresponding 
strong formulation is 

d t u = Au in £2 

pd t u = P Apu + f{u) — d v u onf, 


where the normal derivative d v u figures as the Lagrange multiplier corresponding to the constraint 
Y = yu. 
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2. 3.2 Allen-Cahn equat i on wit h dynamic boundary conditio ns. The following problem is studied 
in Gal & Grasselli ( 2008 ): Liero ( 2013 ): Colli & Fukaol d2014l) and further references therein. Given 
potentials W, Wp : K. —► K. such as a double-well potential W(u) = ( u 2 — l) 2 , the Allen-Cahn equation 

d t u=Au — W'(u) in £2 

is considered with a dynamic boundary condition like in (12.15b . 


fi d,u = j3 Arif — Wp(u) — d v 


on r. 


In Lierol ( 2013 1 this dynamic boundary condition is derived as a scaling limit in a vanishing boundary 
layer approximation. This problem fits into the framework (12.131 ) with (•, •) and a(-, •) as in the previous 
example. For the energy functional 


£(>') = il|Vv|| 2 2(i2) + ^||V r vi| 2 2(r) +lT(v) + Wr(7v), 


ve v, 


(2.16) 


where V = {v£ // 1 (£2) : yv G H 1 (F)} as before, the weak formulation can be viewed as the Li (22) ® 
L 2 (r) gradient flow with respect to the weighted inner product (12.41) : 

(d,u,v) = — E\u)v VvGV. 


2.3.3 Cahn-Hilliard equation with dynamic boundary conditions. Given potentials W.Wp : ffi. —> R, 
the Cahn-Hilliard equation 

d t u = Aw 

, in £2 

w = Au — W{u) 

is considered with the boundary conditions 


/I d t u 
d v w 


Apu — Wp(u) — d v u 

0 


on r 


in 


Cherfils et al. ( 2010l) : Kenzler et al ( 2001 ); Gall ( 2008 ): Racke & Zh eng (2003) as a prototype model 


for the influence of the boundaries on the process of phase separation. 

The corresponding weak formulation is the H~ l (£2) ©/^(.T)-gradient flow of the energy functional 
(12.16b . with the inner product 


(m,v) = (u,v) H -l {a) + il{u,v) Ll( p), 


where = (—with A~ 1 v denoting the solution with zero mean of the Poisson 

equation with homogeneous Neumann boundary conditions for the inhomogeneity v' on £2. This weak 
formulation is again of the form (12.13b with a(-, •) as in (12.14b . 

To our knowledge, this problem is the only parabolic problem wit h dyna mic boundary conditions 
for which a numerical analysis has been carried out: in Cherfils et al. (120id) the Elliott-French finite 
element space discretisation and the implicit Euler time discretisation are studied. 

The above model can be viewed as coupling Cahn-Hilliard in the bulk and Allen-Cahn on the sur¬ 
face. The Cahn-Hilliard / Cahn-Hilliard coupling corresponding to H = H ^(£2) ('\)H '(T) is equally 
of interest; cf. Goldstein et al (2011). 
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3. Spatial semi-discretisation on polygonal domains 

We study the finite element semi-discretisation in space of parabolic equations with dynamic boundary 
conditions in the framework of the previous section, first for linear problems with constant coefficients, 
then for linear problems with space- and time-dependent coefficients, and finally for semi-linear prob¬ 
lems. Mass lumping is also studied. In this section we assume that the domain £2 is polygonal, so that 
no boundary approximation is required. 

3.1 Spatial semi-discretisation of linear problems with constant coefficients 

3.1.1 Galerkin semi-discretisation. For a finite dimensional approximation space V), C V, the Galerkin 
semi-discretisation of the abstract parabolic problem ( 12.11 ) determines u /, : [0, T) —> V), such that 

(u h (t),v h )+a(u h (t),v h ) = (f(t),v h ) Vv/jSV/, (0 <t^T) 

(«/, (0), v h ) = (m 0 , v h ) V v h £ V h . 

Representing u/ft ) with respect to a basis <f»i, (fh,..., (p^ of V* as 

N 

M/,( 0 = 

(=1 

the time-dependent coefficient vector u(f) = (u\(t),U 2 (t),...,UN(t)) T then satisfies the system oflinear 
ordinary differential equations 

Mu(t) + Au(f) = b(f), 

with the symmetric positive definite mass matrix M and the symmetric stiffness matrix A having the 
entries 

m ij = (<Pj,<Pi), a ij = a((p j ,(p i ), (i,j=l,2,...,N), (3.2) 

and with the load vector b(f ) with entries 

bi(t) = ( f(t),(pi ) (i = 1,2,... ,N). 

There is the semi-discrete energy estimate of the same type as dm obtained by the same proof, 

l«*(0| 2 + ||m*(j)|| 2 ds < c 2 " (|M*(0)| 2 + ^^ llrwil^dj) - (3.3) 

3.1.2 Linear finite elements on polygonal domains: first-order error bounds. We discuss in detail 
the case of linear finite elements for the heat equation with dynamic boundary conditions ( 12.51 ). the 
extension to higher-order finite elements being straightforward. Using the Ritz projection corresponding 
to the bilinear form a(-, ■) given by (12.31 ) we show optimal-order error bounds for the spatially discrete 
solution. For the ease of presentation we assume that K > 0 in (12.31 ) throughout this section, so that 
«(•,•) is an inner product on the closed subspace V of H l (Q). We then consider the norm ||v|| 2 = a(v, v) 
on V. The inner product (•, •) on H = L2(£2) ©/^(F) is given by (12.41) . It induces the norm |v| 2 = (v, v) 
on H. 

We consider a family of quasi-uniform triangulations of the domain £2 parametrised by the maximal 
meshwidth h. The corresponding finite element space Vj, C V is spanned by continuous, piecewise 
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linear nodal basis functions <pi, (pi ,..., (Pn that are continuous on £2 and linear on each finite element 
and, for each node Xk, satisfy (pj(xk ) = V We note here that the restrictions of the basis functions to 
the boundary form a basis over the boundary elements. 

The basic tool for proving error bounds is the Ritz projection R/,: V -A V/, with respect to the elliptic 
bilinear form «(•,•) of (12.3b . which is the n-orthogonal projection defined by 

a(R h u,v h ) = a(u,v h ) Vv* € V h . (3.4) 

Lemma 3.1 The error of the Ritz projection for the elliptic bilinear form (12.3b satisfies a first-order 
bound in the energy norm, 

\\V(u-R h u)\\l 2(n) +p\\Vr(u-R h u)\\l 2(n + K\\y{u -R h u)\\l 2(r) 

< [\\u\\~ H 2^ +/3||7 M llH 2 (r))’ 

where the constant C is independent of h and u e 11 2 (£2) with \/~pyu G H 2 (r). 

Proof. With ||v|| 2 = a(v,v) we have for the error of the Ritz projection 


\\u-R h u\\ = min \\u-v h \\ < \\u-I h u\\, 

VfcGVfc 


where /;, denotes the piecewise linear finite element interpolation operator. Using the standard interpo¬ 
lation estimates in the polygonal domain £2 and on its boundary F, we obtain 


\\u-I h u \\ 2 = \\V(u-I h u)\\ 2 ^ w +p\\V r (u-hu)\\ 2 ^ {r) + K\\y(u-I h u)\\ 2 ^ {r) 

< II V(m - Ihu)\\ 2 Ll[a] + P || V r (M - hu)\\l 2{r) + Kc 2 n || u -hu\\ 2 H i {n) (3.5) 

< (i + Kc 2 n )Ch 2 \\u\\ 2 H2{n) + pCh 2 \\u\\ 2 H2{rr 

which yields the stated result. □ 

We note that the order in h increases from 1 to m if finite elements of degree m are used. 

Using the above result for the Ritz map we prove a first-order error estimate in the natural norms for 
the spatial semi-discretisation of (12.51) . 

Theorem 3.1 If the solution of the parabolic problem with dynamic boundary conditions (12.5b is 
sufficiently regular, then the error of the semi-discretisation (13.1b satisfies the first-order bound 


I MO u{t)\\ L2(n) +^i\\y{u h {t) w(0)llz, 2 (r) 

+ J q (j| V M0 - w(0) Hl 2 (C!) + P II' v r MO - «(0) II l 2i r) + K ll 7(«/<(0 - «(0)llL ( r)) ds < c/i 2 

(3.6) 


for 0 < t ^ T. where the constant C is independent of h. but depends on T. 


Proof. Using the error bound of the Ritz projection given in Lemma |3~T1 the proof uses the stan¬ 
dard argument based on the energy estimate; cf. Thomee ( 2006h . We include the short proof for the 
convenience of the reader and as a later reference point. We use the decomposition 


m/j — u = ( m /, — Ri,u) + ( Ri,u — u) 
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where the second term can be estimated using Lemma [XT] Let us denote the first term in the error as 
e/,(t) := H/,(f) —Ri,u(t). Using the definition of the Ritz projection gives us the error equation 

(eh{t),v h ) + a(e h (t),v h ) = (u(t)-R h u(t),v h ) Vv A € V/,, 

and the energy estimate (13.3b (with a = 1 and c = 0) yields 

\e h (t)\ 2 + f\\e h (s)\\ 2 ds^\e h (0)\ 2 + f \\u(s) -R h u(s)\\lds. (3.7) 

Jo Jo 

Since m/,(0) = P/,u( 0) is the //-orthogonal projection of u( 0), we have 

MO)I ^ |P/,m(0)-m(0)| + |m(0)-R/,m(0)| <2|m(0)-R/,m(0)|. 

Moreover, both the dual norm || • ||* and the //-norm | • | are weaker than the V-norm || • ||: C _1 ||v||* ^ 
|v| SC C||v|| for all v' G V. By Lemma IXT1 both terms on the right-hand side of the energy estimate are 
thus 0(}i 2 ) if the solution u is sufficiently regular. □ 

Remark 3.1 The Ritz projection can be defined for the general case KG R. By introducing the positive 
definite forma*(-, •) := a(-, •) +c (•, •) (c is from the Garding inequality), and defining the Ritz projection 
with respect to this modified form. 


a*(R h u,v h ) = a*(u,v h ) Vv* e V h , 

error estimates for the Ritz map and the semi-discrete error bounds can be shown as above. 

3.1.3 Second-order error bound in the case ofWentzell boundary conditions. We consider the spatial 
semi-discretisation of problem (12.5b with [5 = 0 on a polygonal domain using linear finite elements. We 
assume again k > 0 for ease of presentation. We will show second-order error bounds in the L? (£2 ) CD 
// 1/2 (r) norm, using an unusual variant of the Aubin-Nitsche duality argument. We need the following 
H 2 -regularity condition. 

Condition 3.1 There exists a constant C 2 <°° such that for every G LifTl) and gp G // 1 ,/2 (/~), the 
weak solution (j> of the Poisson equation with Robin boundary condition 

—A<t> = go in Q 

^ (3.8) 

d v (j) + K(j) = gp onf 


is in H 2 (Q ) and is bounded by 


ll h 2 ( q ) 


€ 


Ci (JM11^) + Hgrlltf i/ 2(r) ) • 


(3.9) 


Remark 3.2 Condition 13. II is known to be satisfied in the case of a smooth domain see Taylor 
(201 lj). We would expect that it also holds for convex polygonal domains, but we are not aware of a 
reference for such a result. 


The following estimate then holds for the error of the Ritz projection. 
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Lemma 3.2 If Condition l3.1l is satisfied, then the error of the Ritz projection (13.4b corresponding to the 
bilinear form (12.3b with /3 = 0 satisfies the second-order bound 

\\u-Rhu\\l 2( n)+l^\\y( u -Rhu)\\l-i/2 {n ^Ch 4 \\u\\l 2(n) , 

where the constant C is independent of h and u £ 11 1 (Pi). 

Proof. In the spirit of the Aubin-Nitsche duality argument we consider the elliptic problem, for 
g = u- R h u, 

—Aty = g in Q 

d v <p + K(j> = p.(-A r +iy l/2 yg onT, 

which has the weak formulation 

«(<M = {gP’)L 2 (Q)+ Vv€H l (Q), 

where 

{yg,yv) H - i/2 (r) = ((-Ar+iy^yg,^)^. 

Choosing v = g = u — If ,it we obtain 

llsllL 2 (£!)+Ml|y,?lltf-i/2 (r) =a(<l>,g) =a(<l>-R h <l>,g), 

using the Galerkin orthogonality (13.4b in the last equality. This expression is further estimated, using 
subsequently the Cauchy-Schwarz inequality. Lemma IXT1 Condition 13. II and the relation between the 
H l / 2 (r) and H~ l l 2 (r) norms: 

a(<t> —Rhty,g) = a($ — Rhty ,u — Rhu) 

\\<l> -Rh<l>\\ ■ \\u -Rhu\\ 

^ C^ll^ll h 2 (q) ' Ch\\ Ll \\H 2 (n) 

< C'/i 2 ^||g||^ 2 ( I2 ) +/i||(— Ap+I) 1 / 2 yg||^i/ 2(r) ) \\ u \\h 2 (£2) 

= c , /t 2 (||g||^ (i3) +M||yg||^i /2(r) ) 1/ lM|| H 2( i3 ), 

and dividing through yields the stated result. □ 

Lemma l32l yields the following error bound for the spatial semi-discretisation of ( 12.5b with /3 = 0. 

THEOREM 3.2 If the solution of the parabolic problem with dynamic boundary conditions (12.5b with 
/3 = 0 is sufficiently regular and if Condition 13. II is satisfied, then the error of the semi-discretisation 
(13.1b with starting value m/,(0) = R h u( 0) satisfies the second-order error bound 

||M ft (f)-M(0||^( i3 ) + /t||y(n*(f)-M(f))||2_ 1 / 2(r ) <C/i 4 (3.10) 

forO ■Sj t -Sj T , where the constant C is independent of h. 

Proof. We return to the proof of Theorem l3. 1 l and bound Ri,u(t ) — u(t ) using Lemma l32l In the energy 
bound of the error (13.7b we note that e/,(0) = m/,(0) — R/,u(0) = 0 by assumption, and 

ll«A(0ll^(fl)+Mllye*(0ll^-i/2(r) < Il«*(0ll2 2 (fl) + ^l|y«fc(0lli 2 (r) = MOI 2 - 
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Since 

(vv,v) = (w,v)i2( fl ) +ju(yvv,yv) L 2 (r) 

^ 11^11^2(^2) + M||7 vv ll//-i/2(r) lly' , llfl 1 / 2 (r) 

^ c {\\ w \\ l 2 ( q ) + A t ll7 w ll^-i/2( r )) ll v llffi(fl)> 

we obtain 

IMI* ^ c (\H\l 2 (n) + v\\yM\ 2 H -i/2 in ) ', 

which for vv = u{s) — R^ut^s), for 0 ^ s ^ t < T, is bounded by 0{h 2 ) by Lemma [Ol Hence (13.71) 
together with the above estimates yields the result. □ 


3.1.4 Second-order error bound in the case of bulk-surface diffusion coupling. In the case of > 0 
in (12.51 ) we use the following H 2 -regularity condition. 

Condition 3.2 There exists a constant Ci < °° such that for every go € Lt(X 2) and gp c- H l ! 2 (T), the 
weak solution 0 of the elliptic problem, with /3 > 0 and K > 0, 


~A<j)=gi 2 in PI 

<?V0 - Mr0 + K<j) = gr on T. 

is in H 2 (Q), has trace in H 2 (T) and is bounded by 

Il0llff 2 (t2) + H^ll// 2 (r) ^ C 2 (\\g n \\l 2(n) + ll^r|lz 2 (r))- 


(3.11) 


(3.12) 


Remark 3.3 It can be shown that this condition is satisfied for smooth domains PI, adapting the 
standard H 2 -regul arity pr oof for the heat equation with Dirichlet or Neumann boundary conditions 
as given, e.g., in ( Brezisi. 2011 . Section 9.6). We are not aware of a suitable reference for (convex) 
polygonal domains. 


Lemma 3.3 If Condi tion l3.2l is satisfied, then the error of the Ritz projection ( 13.41 ) corresponding to the 
bilinear form ( 12.31 ) with > 0 satisfies the second-order bound 


h- R iA\l 2 (n) + ^\W-Rhu)\\ 2 L2(r) < Ch 4 (IMI h 2 (q) + IMI^ 2 (r))’ 


where the constant C is independent of h and u £ H 2 (Pi) with yu £ fi 2 (r ). 


Proof The proof uses the standard Aubin-Nitsche duality argument with go = u Rh u and gr = 
y(u — Rhu) in (13.1 II ): cf. the proof of Lemma [X2l We omit the details. □ 

Together with the proof of Theorem l3.ll Lemma l33l vields the following error bound for the spatial 
semi-discretisation of (12.5b with /3 > 0. 

THEOREM 3.3 If the solution of the parabolic problem with dynamic boundary conditions (12.51) with 
ft > 0 is sufficiently regular and if Condition 13.21 is satisfied, then the error of the semi-discretisation 
(13.11) satisfies the second-order error bound 

I MO - “( f )lli 2 (i2) + M||7 (m/.( 0 - «W)lli 2 (r) < Ch4 

forO ^ t ^ T , where the constant C is independent of h. 


(3.13) 






PARABOLIC PROBLEMS WITH DYNAMIC BOUNDARY CONDITIONS 


13 


3.1.5 Mass lumping. In many situations it is convenient to replace the mass matrix M of (13.21) by a 
diagonal matrix. This can be achieved by replacing the L 2 (Q)®L 2 (r) inner product (•, •) of (12.41) by a 
suitable quadrature. 

We follow the framework of lumped-mass methods as presented in dThome 5 120061 Section 15). We 
use the analogue of the trapezoidal rule in d and d— 1 dimensions: 


j d +l 2 ^ 

Qe (/) = vol,/ {E ) —— £ f(xj ) and Q e (/) = vol rf _ i (e) - £ f(xj ), 

“ + 1 j= i a j= i 

where vol c / denotes the rf-dimensional volume, E is an element of the quasi-uniform triangulation 
with vertices Xj, and e G r):% is a boundary element. Then we use the above quadratures to define the 
lumped mass approximation of the inner product (•,•) on // = LjiQ) ©/^(L): 

(«, v)lm := Y Qe{uv) + p Y Qe(uv), (3.14) 

Ee3r h eed3r h 


We denote the induced (lumped mass) norm by | • 


The lumped-mass semi-discrete problem determines «/, : [0, T] —> V), such that 

(uh(t),Vh)iM + a(uh(t),Vh) = (f(t),Vh)iM Vv/i £ Vh (0 <t^T) ( 315 ) 

(m/,(0),v/ ! )lm = («o,p/i)lm Vv/j G Vh 

Then the nodal vector u(f) satisfies the corresponding linear ODE system 

Mu(l) + Au(f) = b(f), (3.16) 

where the new mass matrix, which we denote again by M, has the entries 


m ij = ((Pj,<Pi) lm, and bi(t) = (f(t),(pi)uA 


The new mass matrix is diagonal, since the product (pj(pj vanishes at all nodes with j / i. The stiffness 
matrix A is defined as before, see (13.21) . The initial values are now chosen such that ui,(xj, 0) = u(xj, 0) 
in each node of the triangulation. 

Before proving error estimates for the semi-discretisation with mass lumping we formulate an im¬ 
portant technical lemma. Here || • || is again the norm induced by the bilinear form a(-, •) of (12.3k 

Lemma 3.4 If /3 > 0 in (12.41 ). then the quadrature error S{vh-Wh) = (v*,vv*)lm — (v/,,w/,) is bounded 
by 

K(v/„w/,)| < Ch 2 ||v/,|| ||w/,|| for v h ,Wh G V h . 


Proof We separate the mass lumpin 
proof of Lemma 15.1 in Thomee (200( 


errors in the bulk and on the surface as £ = £q + p £‘r ■ The 
shows that 


Wn{vh,Wh)\ ^ Ch ||v7i||//i(r2) Il vv /ill//i(i2) 
\^r(yh,Wh)\ f C/ ! “|| l 7i||//i(r) ll^ll/rqr) 


for v/j, w h G V h . 


(3.17) 


This yields the stated bound for £ for the norm || • || induced by the bilinear form a (•, •) of (12.3k provided 
that p > 0 . □ 
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Theorem 3.4 If the solution of the parabolic problem with dynamic boundary conditions (12.51 ) is 
sufficiently regular, then the error of the semi-discretisation with mass lumping (13. 151) satisfies a first- 
order error bound as in Theorem 13.II If (3 > 0 and Condition 13.21 is satisfied, then there is also a 
second-order error bound as in Theorem l3.3l 

Proof. The proof differs from that of Theorem l3.1l in that the error equation for the lumped mass case 
for eh{t) = u/,(t) —Ri,u(t) contains extra quadrature error terms: for all v/, £ V/,, 

{e h (t),v h ) LM +a(e h (t), v h ) = (u(t) - R h u(t),v h ) - £(R h u(t),v h ) + £{I h f(t),v h ) + ( 4/(0 - f(t),v h ). 

We first consider the case (3 > 0. The energy estimate obtained by testing with v/, = £/,(/) and using 
Lemma 13741 to bound, for w/, = Ri,u(t) — Ihf(t), 

£{wh,e h ) <C/i 2 ||w /l || \\e h \\ < ;j||e/ I || 2 + C7i 4 ||w, ! || 2 

then yield the result with | • |lm in place of | -1. Together with the inverse estimate (valid for quasi¬ 
uniform triangularisations) 

INXC/tVaI for Vh £ 14, 

Lemma [3~4l further provides the equivalence of the norms | • | and | • |lm on V/„ uniformly in h. 

In the case f3 = 0 we use (13. 171 ) and the inverse estimate to obtain (with different constants that are 
all denoted C) 


<?(wh,eh) ^ C/z 2 (||w/ ! || // i( i 2 )||e/,|| // i(j 2 ) + || w/, ||#i (r)ll e /, ll//i (r)) 

^ C/z 2 (||vv/ ! || ff i( fl )||e/ ! || ff i( i3 ) + ||w/,|| ff i( r )C/i 1 lk/)||z,2(r)) 

< Ch{\\w h \\ H i( a) + ||wA|| ff i (r) ) Ik/,11 

< i ll e /il | 2 + C/z 2 (||w' /i ||^ 1{i2) + ll w /ill^i( r) )- 


This is sufficient to obtain the first-order error bound with | • |lm in place of | • |. Estimating the quadrature 
error in the above way still yields the /z-uniform equivalence of the norms | • | and | • |lm ° n V/,. □ 

Remark 3.4 A second-order error bound is not obtained for mass-lumping in the case of Wentzell 
boundary conditions (j3 =0). A second-order error bound does hold, however, if mass-lumping is done 
only in the bulk, but not on the surface. For such a partial mass-lumping Lemma [3741 remains valid for 
13=0. 


3.2 Spatial semi-discretisation of lion-autonomous linear problems 

The results of the previous subsection generalize to the heat equation with non-autonomous dynamic 
boundary conditions (12.101 ) on a polygonal domain. Using the time-dependent Ritz projection corre¬ 
sponding to the bilinear form given by ( 12.1 lb . we show optimal-order error bounds for the 

space discretisation. Since the proofs of these results are very similar to the ones presented above, we 
only focus on the major differences. Again for the ease of presentation we assume K(x,t) > 0 and 
bounded by K, while f3 is bounded by f3, hence the form a(t;-,-) induces the norm ||v|| ; = (/(Lv.v ) 1 / 2 
on V = {v G H l (£2) : Vf-(yi') £ Lt(F )}. The inner product on H= Li(^) © L 2 CO I s given by 

(12.12b . and induces the norm |w| ( = m[t\ v,v ) 1 / 2 on H = L 2 (£l)®L 2 {T). Note that both norms are time 
dependent, but in each family of norms the members are equivalent to each other uniformly in t. 
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We work with the same family of quasi-uniform triangulations as in the previous section, and 
also use the same piecewise linear finite element basis functions. The semi-discretisation of the non- 
autonomous problem is 

m(f,u h (t),v h )+a(t,u h (t),v h )=m(f,f(t),v h ) Vv A eV* (0 <t^T) 

(3.1o) 

m(0;u h ( 0),v A ) = m(0:u Q ,v h ) Vv/, £ V h . 

This now yields a non-autonomous system of linear ordinary differential equations for the coefficient 
vector u(f), 

M(f)u(f) + A(f)u(f) = b(f), 

where M(f) and A (t) are the time-dependent mass and stiffness matrix, respectively. Their entries are 
given as 


mij(t)=m(t;<Pj,<Pi ) and a, 7 (f) = a{t;<pj,<Pi) (i,j=l,2,...,N), (3.19) 

and b(r) is the load vector, with entries 

bi(t)=m(t,f(t),<Pi) (i= 1,2,..., N). 

In order to prove optimal order error estimates we begin by studying the Ritz projection R/,(f) : 
V —> Vh with respect to the elliptic bilinear form a(t ; •, •) of (12.111) . for 0 ^ t ^ T. For every u £ V, the 
projected function Rh(t)u £ V/, is defined as the unique finite element function in V/, such that 

a{t\R h (t)u,v h ) = a(t;u,v h ) Vv h &V h . (3.20) 

The error bound of Lemma [3711 extends to the time-dependent case using the same proof and observing 
the uniformity of the bounds in t. 

Lemma 3.5 The error of the Ritz projection corresponding to the bilinear form (12.1 11 1 satisfies a first- 
order error bound in the time-dependent energy norm for 0 ^ f ^ T, 

\\u - Rh(t)u\\j < Ch 2 ^\\u\\ 2 H 2 {ii) +My u \\ 2 H 2 {r) ), 

where C is independent of h and t. Here, ji is an upper bound of j3 (•, t). 

We have Rj,(t)u(t )) 7 -Rh(t)u(t) in general, but we note the following bound. 

Lemma 3.6 The error of the time derivative of the Ritz projection corresponding to the bilinear form 
( 12.1 It satisfies a first-order error bound in the time-dependent energy norm for continuously differen¬ 
tiable u : [0, T] —> H 2 (Q) and 0 ^ ^ T, 

ll«( 0 -|(^( 0 «( 0 )ll?< c/z 2 (|| M (0ll^2 (i 2 ) +y3||r«(0ll^ co + l|M( 0 ll^ ( ^ ) +j 8 |lr«( 0 ll^ 2 (r) ), 

where C is independent of h and t. Here, /j is again an upper bound of j3(-,f). 

Proof. Differentiating (13.20b with respect to t yields 


a(f,u{t) - j- t (R h (t)u{t)),v h ) = -d,a(t\u(t)-R h (t)u(t),v h ) Mv h £ V h . 
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We then have, omitting the argument t, 

||« — ^{ R hu)\\j = a(t',u— j- t {R h u),u-R h u)+a(r,u- j- t (R h u),R h u- f t (R h u)) 

= a(t\u- j- t (R h u),u-R h u)-d t a(t\u-R h u,R h u- f t {Rhu)) 

= — jj(Ri,u),u —Ri,u) — d t a(t\u — R/,u,RhU — it) — d t a(t\u — Ri,u,u— R/,u )). 

With (12.7b and Lemma l-Ol this is estimated as 

\W-i(Rhu)\\f < ll«-|-(^/ i «)llfV / C/!(||h||^ 2(i2) +^||7H||^ 2(r) ) 1 

+ M[Ch 2 (\\u\\~ H 2 ^ + PWyuWhh^ (ll“ll^ (i2 ) +/ 3 ||y«llH 2 (r)) 

+ M[Vch(\\u\\ 1 H 2^ n s j +^||r«ll^(r)) ll“ — ?[i( R hu)\\t, 

which yields the result. □ 

We then have the following non-autonomous version of Theorem l3.ll 

Theorem 3.5 If the solution of the linear non-autonomous parabolic equation with dynamic boundary 
condition (12. 10b is sufficiently regular, then the error of the semi-discretisation (13.18b satisfies the first- 
order error bound 

IKO) - u (t)\\l 2 (n) + J (jUh{t) - 7«0)) da + J ||' V(u h (s) - u(s)) ||^ (fl) ds 
+ J J j3(.,s)(Vj-n/ ! (s) — V r u(s))~dods + J J K(.,s)(yuh(s) — yu(s)) 2 d(rds ^Ch 2 , 

for 0 < t ^ T, where the constant C is independent of h and t, but depends on T. 

Proof. The error equation for e/, = u /, — R/,u is now 

m{f,e h (t),v h ) + a(t-,e h (t),v h ) = m(t;u(t) - jj(R h (t)u(t)),v h ) Vv h eV h . 

The energy estimate ( 12.9b (with a = 1 and c = 0) thus gives us 

\ eh (t)\f+J o \\eh(s)\\*ds ^ e 2 ^ (\e h (0)\l + ||n(s) - £ s (R h (s)u(s))\\l tS ds > ), 

and the result follows with Lem mas 1.3. 5 l and [Thl □ 

Under appropriate // 2 -regularity conditions, second-order error bounds are obtained in the same way 
as in Subsections l3.1.3l and l3.1.4l 

Mass lumping can be discussed as in Subsection l3.1.5l In particular, the second-order error bound 
of Theorem l3.5l remains valid for mass lumping both in the bulk and on the surface for strictly positive 
/j (x.t), and for mass lumping in the bulk, but not on the surface for the case p = 0. 


3.3 Spatial semi-discretisation of semi-linear problems 


The above techniques readily extend to semi-linear parabolic problems with dynamic boundary con¬ 
ditions. Let us consider the following semi-linear problem, which includes the first two examples in 
Sectionl231 


d,u = An + fn(u) 


p d t u = jiAru — ku + pfr{u) — d v u 


in PI 
in r. 


(3.21) 
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with /I > 0 and ^ 0 , and where the nonlinearities fa : R —> R and fr : K —> R are continuously 
differentiable functions, and hence locally Lipschitz continuous. For u G C(£2) we use the notation 

/(«) = {fn{u),fr(u)) T G C(f2)©C(r) with fa(u)(x) = f a (u{x)) forxeQ and fr(u)(x) = f r (u{x)) 
forx G r. 

The bilinear forms corresponding to the problem are the same as in ( I2.3b -( 12.4b . again taken with a 
positive K to simplify the presentation (otherwise the linear term can be absorbed in fr(u)). The finite 
element semi-discretisation of the semi-linear problem is 

(uh(t),v h ) + a{u h (t),v h ) = (f(u h {t)),v h ) Vv h eV h (0 <t^T). (3.22) 

For the coefficient vector u(f) this yields the nonlinear system of ordinary differential equations 

Mu(f) + Au(f) = f(u(f)), 

where M is the mass matrix, A is the stiffness matrix, both defined in (13.2b . and f(u(f)) is the vector 
with entries /;(u(t)) = (/(«/,(f)), (pt) ( i = 1,2We give the semi-linear analogue of Theorems 

IXTtiPl 

Theorem 3.6 Consider the semi-linear equation with dynamic boundary condition (13.2 lb over a polyg¬ 
onal domain Q C R 4 with d f 3. The following error bounds hold if the solution of ( 13.211 ) is sufficiently 
regular, if Conditions 13. II and 13.21 are satisfied and if the initial data satisfy ||m/,(0) —Rhu( 0)11^^ + 
/i||y(w/,(0) — /?/ 1 n(0))||^ r -j ^ Co/i 4 . Then, the error of the semi-discretisation (13.22b by linear finite 
elements on quasi-uniform triangulations satisfies the first-order error bound 

\\u h (t) +MllrM 0 -M( 0 )lli 2 ( r) 

+ J o (j| v Mr) - «( s )) ||^( fl ) + P || v r(M/.(s) - m(s)) II l 2in + k||t(ma(s) - w(5))Hi 2 (r ) Jds < Cfi 2 , 

(3.23) 

for 0 < t ^ T and 0 <h f /io with a sufficiently small ho. Moreover, there is the second-order error 
bound 

\Wh{t) - u(t)\\l 2in) + l^\\r(uh(t) ~ C 2 h 4 , (3.24) 

where 0 = 1/2 if /3 = 0, and 0=0 if jj > 0 in (13.211) . The constants Ci and Of are independent of h 
and t, but depend on T. 

Proof. The error equation for e/,(f) = «/,(f) — Ri,u(t) reads, with the elliptic bilinear form a of (12.31 ) and 
the inner product (12.4b . 

(eh(t),Vh)+a{e h {t),v h ) = (u(t) -R h u(t),v h ) + (f{u h (t)) - f(u(t)),v h ) Vv A e V h . 

We use again the energy estimate, which now becomes 

\ e h{t)\ 2 + [ \\e h (s)\\ 2 dsif\e h (0)\ 2 +[ \\u(s)-R h u(s)\\lds+ f \\f(u h (s))-f(u(s))\\lds. 

Jo Jo Jo 

As long as ||«( 0 llz., x ,(t 2 ) ^ r an d H M /;( f )llL M (t 2 ) ^ r - the local Lipschitz continuity of f n and fr yields 
(omitting the argument t ) 


II/(«/,) - /(«)IU < L{r)\u h -u | < L{r)\e h | + L{r)\R h u -u\. 
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The maximum-norm boundedness of u /, is ensured by an inverse estimate: 

Il“/i||z«„(i2) ^ \\ u h — IlM\uo{S2) + WhP\\u,{n) 

< Ch d /~\\uh — huW^a) + IMIm* 2 ) 

< Ch d/ -\\u h -u\\ L2{£i )+Ch d /~\\u — h^W^Q) + 

< CC 2 h 2 - d l 2 +Ch 2 - d l 2 \\u\\ H 2 (a) + ||H|| M fl) 

^ 1 + IMII»U2) 

for sufficiently small h, where the last but one inequality holds true as long as (13.24b is valid. The error 
bounds for the Ritz projection given by Lemmas |3.1113. 3 1 and a Gronwall inequality conclude the proof. 
□ 

Mass lumping can be studied as before, yielding the semi-linear analogue of Theorem l3.4l 


4. The Ritz map for non-polygonal domains 

For smooth domains the polygonal approximatio n of the bulk and the surface requires extra care in the 
error a nalysis of spatial discretisations; cf., e.g., Dziuk ( 1988 ): Dziuk & Elliotd ( 20071) : Elliott & Rannerl 
(20131). In this section we will discuss spatial semi-discretisations of problems with dynamic boundary 
conditions on a smooth domain £2 with boundary F. 


4.1 Preparation: Lifts and their approximation estimates 


In this subsection we d escribe the setting and recall some approximation results from Dziuk (1988); 
Dziuk & Elliott ( 2007 ): Elliott & Ranneii ( 20131) . 

The smooth domain £2 is approximated by a polygonal domain £2/, with boundary surface I), := d£2 /„ 
triangulated with a mesh size h that is assumed sufficiently small in all the following. In particular, we 
require that for every point i£l/, there is a unique point p G F such that x — p is orthogonal to the 
tangent space T p r of F at p (see Figure [!}• We assume that the vertices of I/, are on F, i.e., T/, is 
an interpolation of F. We consider a family of quasi-uniform triangulations parametrised by the 
maximal meshwidth h. 

Following Dziuk (1988), we define a lift of functions v/,: T/, —► R to v l h : F —► R by setting v l h {p) = 
v/j(x) for p G F, where x € F, is the unique point on F/, with x — p orthogonal to the tangent space 


T p r. As in lElliott & Rannei ( 2013 ). we define a lift of functions Vh : £2h —> R to v[ : Q — > K by setting 


v l h (p) = v h {x) if x G £2/, and p G £2 are related as in Figure [Q see Elliott & Ranner (2013) for a precise 


'h 

formal definition. Note that both definitions coincide on F. 


The finite element space V), ^ H l (£2) corresponding to :% is spanned by continuous, piecewise 
linear nodal basis functions on £2/,, as in Section 13.1.21 We note here that the restrictions of the basis 
functions to the boundary I/, again form a basis over the approximate boundary elements. We denote by 
Vj t = {v l h : v;, G 14} C F the lifted finite element space. 

Lemma 4.1 For v G H 2 (£ 2), such that yv G // 2 (F), we denote by //, v G V l h the lift of the finite element 
interpolation f t v G V),. Then the following estimates hold: 

(i) Interpolation error in the bulk; see [Bernar d3 dl989b : lEmott & RanneJd2013b : 


l' , - / Av|| i 2 (fl) +/l||V(v-4v)|| L 2 (fl) <C/l 2 ||v|| fl 2 (fl) . 
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Fig. 1: (a) shows the lift for surface functions: v l h {p) = v/, (x) ; (b) shows the lift for bulk functions: 

v l h (p) = n(x). 


(ii) Interpolation error on the surface; see 


iDziukl rtl988h : 


ll7(v-4v)|| L 2(r)+ / t||V r (v-4v)|| L 2 ( r) < C7r 2 |Nlr4(r)- 


We introduce the discrete counterparts of the bilinear forms (12.31) and (12.4b . respectively: 

a h {u h ,Vh)= / Vi^-Vv/jdx+K- / {Yi,u h )(Y h v h )dG h + p / V r v h da h (re > 0, y3 ^ 0), 

Jn,, Jn Jn 


m h {uh,Vh)= / u h v h dx + p {YhU h ){YhVh)&Oh 
J Clu J Tu 


i n h 

where Yh denotes the trace operator onto //,, and V j- is the discrete tangential gradient. 


(4.1) 

(M > 0), (4.2) 


Using the lift and its properties (see Section 5 inlDz iuk & Elliott (2013) and Section 4.3 in Elliott & Ranner 
( 20131) for the surface and the bulk, respectively), the following estimate of the geometric errors is ob¬ 
tained. 


Lemma 4.2 For any v/,, w/, £ V/ ; we have the estimates 
\a(v l h , w h)- a h(vh,w h )\ < C/i||Vv{,|| l2(b , ) \\Vw l h \\ L 2 {B , h) 

+ C/i 2 (j|Vv'|| ffl(n) ||Vw'|| L 2 (i2) +j3||V r vi|| i 2 (r) ||VrwJ,|| L 2 ( r ) + ^IItV/JIl 2 ^) WWhWilp)) > 
\m{^ h ,w l h )-m h {v h ,w h )\ < C/tllv^ll^/) 

+ C/!“(^||v/ ) || i 2 (i2 ) ||w/,|| i 2 ( fl) + 411T 1 '* 11/3(0 llMllo(o)’ 


where B | denotes the layer of lifted elements which have a boundary face. 

This result can be shown in the same way as Lemm a 6.1 in Elliott & Rannei ( 20131) . Lemma l4~2l is 
the extension of Lemma 5.5 in Dziuk & Elliott! ( 20131) to the surface-bulk case. As a consequence of 


this lemma we also have the /z-uniform equivalence of the norms \\v l h 


and 


\Vh\h 
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The following lemma provides a key estimate. 


Lemma 4.3 ( Elliott & Ranner ( 20131) . Lemma 6.3) For all v £ H 1 (£2) the following estimate holds: 


IMIl 2 (b') ^ Ch? IMI# 1 ^). 


(4.3) 


4.2 The Ritz map for problems with time-independent coefficients 

The semi-discretisation of the parabolic problem with dynamic boundary conditions (12.5b . with a non¬ 
poly gonal domain determines u/,: [0, T] —> V* such that 

m h (uh(t),v h ) + a h (u h (t),v h ) = m h (f(t),v h ) Vv* £ V h (0 < t f T) ^ 

m(uh( 0),v/,) = m h (u 0 ,v h ) Vv/j £ V h . 

The Ritz map ff, : V -A V l h is defined by first determining R/,u £ V), for u £ V via 

a h {RhU,Vh) = a(u,v l h ) \fv h £ V h (4.5) 

and then setting Ri,u := ( Rf,u) l £ Vj r 

The Galerkin orthogonality does not hold in this situation, just up to a small defect. The following 
estimate is obtained from Lemma l4~2l 

a(u-R h u,v l h ) = a(u,v l h ) - a(R h u,v l h ) = a h (R h u,v h ) - a{(R h u) 1 ,v l h ) 

^ Ch\\VRi,u\\ L 2^ B ij \\^v l h \\ L 2^ B i^+Ch-^\\VR i,u\\ h i^ ||VvJ,|| L 2 ( fl ) 

+ P\\VrRi,u\\i2(r) ll v r* ; /!llL 2 (r) + K \\Y R h u \\L 2 (r) IIllz, 2 (r)^ • (4-6) 

We prove the following error estimate for the Ritz map. We again set || v|| 2 = a(v, v). 

Lemma 4.4 The error of the Ritz map for the elliptic bilinear form ( 12.31 ) on a smooth domain satisfies 
the first-order error bound in the energy norm 

\\u-R h u\\ 2 < Ch 2 (jM|tf 2(i2) + p ||yw|lH2 (r) ') , 

where the constant C is independent of li f ho (ho sufficiently small) and u £ ll 2 (Q) with \fPyu £ 
H 2 (T). 

Proof Using the bound (14. 6I ). we start by estimating as 

|| u — Ri,u\\ 2 = a(u —Ri,u,u —Ii,u) +a(u — Ri,u,Ii,u — R/,u) 

< \\u-R h u\\ \\u-I h u\\ 

+ Ch\\VR h u\\ L2{B , h) \\V(I h u-Ri,u)\\ L 2 [B ij+Ch 2 \\R h u\\ \\I h u-R h u\\. 

These three terms need to be estimated in a suitable and careful way. We use the interpolation estimates 
from Lemma |4~I1 to obtain for the first term 

||k-R /! h||||w-4w|| ^ ^||m — Ri,u\\ 2 + Ch 2 + P ll7 M llH 2( r ) ^ ■ 
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For the second term we have 

Ch\\WRhu\\ L 2^ g i^ ||V(//,m — Rh u ) 

^ Ch(\\V{u -R h u)\\ L 2 (B i h) + ||V M || l2K) ) (||V(m-/ a m)|| Zi 2(^) + ||V(m — Rhu)\\i2( B r^ 

< C/i(2||V(m-/? a m)||£ 2( ^ ) + ||Vm||^ ( ^ ) + ||V(m-4m)||“ 2(5 , ) ) 

where in the last estimate we used Lemma [431 and the interpolation estimate. 

The last term is estimated, using Young’s inequality and the interpolation results of Lemma l4~Tl as 

Ch 2 \\Ri,u\\ \\Ii,u — R/ t u\\ 

< Ch 2 [\\u-R h u\\ + ||m|Q [\\u-I h u\\ + ||m-R/,m||) 

< Ch 2 (2\\u-R h u\\ 2 + \\u\\ 2 +\\u-I h u\\ 2 ) 

< Ch~{2\\u—Rhu\\~ +C{\\u\\ 2 h i^ + /3||7“llffi(r)) +^^~(\\ u \\h 2 (£ 2 ) + P\\y u \\H 2 (r))')■ 

Reinserting all these into the estimate above, absorbing the terms ||n — Ri,u\\ 2 using an h ^ h o with a 
sufficiently small h o, we obtain the stated result. □ 

The second-order error bounds of the Ritz map can be shown analogously as in Section^ Note that 
Conditions 13. II and 13.21 are satisfied for smooth domains Q.. 

Lemma 4.5 The second-order estimates of Lemmas 13.21 and [3731 hold for smooth domains, for h ^ ho 
with a sufficiently small ho- 

4.3 The Ritz map for problems with time-varying coefficients 

In the time-varying non-polygonal case the Ritz map Rift) : V —> with respect to the bilinear form 

( 12. lib , for 0 t T, is defined via R/,(f )n € 14 defined by 

a h {f,Rh(t)u,v h ) = a(t;u,v‘ h ) \/v h € V hl (4.7) 

and setting Rh(t)u := (R/,(t)n) / € Vj r 

Lemma 4.6 The error bounds for the Ritz map and its time derivative given in Le m m as [3. 5 1 a 1 1 d [3761 also 
hold in the case of non-autonomous problems and a smooth domain £2. for h <C ho with a sufficiently 
small ho- 

Proof. In order to prove the error estimates of the Ritz map in the non-polygonal time-varying case 
we have to make sure that the constants in the estimate of the geometric error a/,(f;M*,v/,) — a(f;w/,,V/,) 
(from Lemma 1431 ) are independent of t. This holds because of the smoothness of the domain £2. and 
since we assumed uniform bounds on the coefficient functions. Therefore the first order error bound is 
obtained with the proof of Lemma l4~4l 

The error of the time derivative of the Ritz map is shown similarly as in the proof of Lemma 13.61 
We start by differentiating the definition of the Ritz map (14.71) . and obtain 

a[t\u(t),v l ^j \a h (t\R h {t)u(t),v^j - d t a{t\u(t) 
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Using (g \\’h) 1 = ^ v l h , we estimate 

d , 2 / d 


dr 


(Rhu) = a(^t;u — — (R/,w), m — Rhii^j +a(t',ii— — — — (R/,m)^ 

^ a(t;u — — (R/,m),m — Rhu'j ~ d t a[t\u — Ri,u,Ri,u -(R/,m)^ 

Qh ( r ' ^^ /,M ) : ^ hl1 ~ <Ir _ ( f: ^ (Rhu),RhU - ^(R/ju)') 
d t a h [t\R h u,R h u- —(Rhu)^ - d t a(t;R h u,R h u - ^(R/,«)) 


The first two terms are estimated analogously in the proof of Lemma 13.61 The other two terms are 
estimated separately, like in the proof of Lemma 14.41 but here we use the error estimates for the Ritz 
map instead of the interpolation estimates. For the third term we have by the time-varying version of 
Lemma 14.21 

a h{t\^-{Rhu),RhU- —{RhU^j -a(t\^-(R h u),R h u- 
^ Ch(2\\u- — (Rhu)\\ t +Ch\\u\\ 2 H 2^ n s ) +Ch~\\u\\~ H 2^ £1 ^j 

+ Ch-(2\\u- — (Ri,u)\\ t +C(||m||^ 2 ( i q) + || +C/i 2 (||ii||^ 2 ( fl ) + II V^7“ll H 2 (r))\ 


For the fourth term we note 


d t a(t;v,w) = J $ V r v • V r vvd(7 + j K(yi')(yw)d<7, 

which only contains boundary terms with bounded coefficient functions. As in Lemma l4~2l we obtain 
the bound 

|^a(f;v^w{,) —^a A (f;v A ,w A )| <C/i 2 ^||v / ^V r v[|| t 2 (r) ||v / ^Vrw*|| i 2 (r) + ||\/K7v{,|| L 2(r) llv^MlIz^r))- 

Therefore we have for the last term, similarly as for the previous one, 

d t a h (t\R lt u,R lt u- ^(R h u)^j - d,a(t\R h u,R h u- 

.. d ,, 

^ Ch \\Ri,u\\ t ||R/,m — — (R/,«)|| f 


€ 


Ch~[\\u— —(Rhu)\\ t + C(\\u\\~ Hl ^ Q) + \\VpY*\\ 2 ti2(r)) + II V^'Hl^r)))- 


Absorbing the terms C/i||m— ^-(R/,m)|| with a sufficiently small h in the left-hand term, the result 
follows. □ 


4.4 Error bounds of the semi-discretisation 

With the obtained error bounds for the Ritz map and the geometric error bounds, all the proofs of error 
bounds given in Section[3]now extend directly to the case of smooth domains. We summarize this in the 
following theorem. 
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THEOREM 4.1 The semi-discrete error bounds of Theorems 13.11 - [3761 also hold for smooth domains, 
for h ^ /zq with a sufficiently small /zq. 


5. Time discretisation 


5.1 Standard implicit integrators 

Since the parabolic problems with dynamic boundary conditions are cast in the same abstract setting of 
parabolic problems in which time integration by backward difference formulae (BDF) or algebraically 
stable implicit Runge-Kutta methods (such as the Radau IIA methods) has been analyzed before, results 
on the s emi-discretisatio n in time by these methods can be applied dire ctly; see, e.g.. lAkrivis & Lubich 
(12015h : Crouzeixl ( 1 975l) : Lubich & Osterma'nnl ( 1995 ): Savare ( 19931) for stability and error analyses 
for linear, semi-linear and quasi-linear problems in an abstract setting that applies also to the problems 
considered here. Together with error bounds of the defects obtained on inserting the Ritz projection 
of the exact solution into the spatial semi-discretisation, such as derived in Sections 3 and 4, one then 
obtains error bounds for the full discretisation from the known results. 

As an illustration of this procedure, we consider the k-step BDF time discretisation / finite element 
space discretisation of the linear non-autonomous problem (12 .10b : with the bilinear forms (12 .111 ) and 


(12. 121) . the finite element space 14, a time stepsize T > 0 and given starting values , 
determine zzj| £ 14 for n ^ k with t n = m < 1 from the equation 


1 


7i>''' ’ M /r 


£ 14, we 


™{t n \d?u n h ,v h ) + a(t n ;u n h ,v h ) =m{t"\f{t n ),v h ) Mv h £ V h , 


(5.1) 


where denotes the k-th order backward difference approximation to the time derivative. 


dy h = l -isjur j 

r t =0 


with the coefficients 8j given by the generating polynomial YJj=o&j& = Tle=\ |(1 — CY- I n matrix- 
vector form, this is the implicit scheme 

+ A(f"))u" = b(f") - M(f")i £ SjU n - j . 

We then have the following fully discrete analogue of Theorem l3.5l 

THEOREM 5.1 If the solution of the linear non-autonomous parabolic equation with dynamic boundary 
condition (12.10b is sufficiently regular and if the starting values are sufficiently accurate, then the error 
of the full discretisation (15.1b with linear finite elements and the k-step BDF method with k ^ 5 satisfies 
the following error bound: 

K-<n\\l i n ) + f M(^ n )(r^-7«(f n )) 2 dc7 + Tt||V( M / J - M (^))||^ 

Jr j=k 

+ t E / i 3 (--f ; )( v rM/ ! -V r zz(f ; ')) 2 d(7 + T^ / <;t j ){yu{ - yu(t j )fdo < C(h + x k ) 2 , 

j=k Jr i=k Jr 


for ti ^ k and t — u'Z ^ T , where the constant C is independent of /i, t snd t , but depends on T. 
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Proof. The error equation for e^ = u )j — Ri,(t n )u(t n ) is 

m(t n ; d, T e n h ,vi,) + a(t n \e n h ,v h ) = m(t n \d’i,,v h ) Vv h G V,„ 


where 

d n h = (u(t n ) - ( dfu){t n )) + d?(u -R h u) (f"). 

Here the first term is bounded, for example in the norm ||v|| 2 = a{t\v.v) or |v| 2 = m(t\v,v), by 0(X k ) 
for a temporally smooth solution u. The second term is a linear combination of k shifted difference 
quotients 


1 


u — Riju) (t n ■*) — (u — R/,u) (t n 




1 d 
dt 


{u-R h u){t n - j ~ l + 0T)d0. 


Using Lemma 13761 we thus obtain 

\\d n h \\ t n^C(T k + h). 

The following stability estimate is proved in lLubich et al. 1 ^20l3b by testing with Vh = e” — tlel , 1 with 


T7 = 0, 0, 0.0836, 0 . 2878 0.8160 for k — 1,2,...,5, respectively, using results of iDahlciuistl (119781) : 
Nevanhnna & Odeh ( 1981 ) and the bounds (12.61) and (12.7b : for t " L T, 


e n h \l- 


+ ^i\\ e X \\dj 


j=k 


i t j+C max 

’ r (KKk- 




Using Lemma [3751 and recalling the definition of the bilinear forms a and m of (12.3b —( 12.4b then yields 
the result. □ 

We note that the order in h increases from 1 to m if finite elements of degree m are used. 

The result for the semi-linear equation, Theorem 13.61 is extended to the BDF full discretisation in 
the same way under the mild stepsize restriction x k f Ch d/ ' 2 , which allows us to show the L°° bound of 
the fully discrete solution in the same way as in the semi-discrete case. 


5.2 Exponential integrators 

Exponential integrators have become attractive for use with parabolic problems in recent years_since they 
allow for large time steps while using just matrix-vector multiplications; see the review by Hochbruck & Ostermannl 
(2Q10). While it is not obvious from the outset how to use exponential integrators for parabolic prob¬ 
lems with dynamic boundary conditions such as (1231 the above weak formulation and finite element 
discretisation with mass lumping yields the system 


Mu(f) + Au(f) = b(f) 

with a diagonal positive definite mass matrix M. With A = M _1 / 2 AM _1//2 and b(f) = M _1 / 2 b(r) and 
y(f) = M | / 2 u(f), the transformed system becomes 

y(t) = -Ay(t) + b(t), 


to which an exponential integrator is readily applied. Here we consider for simplicity of presentation 
just the exponential Euler method 


y" + T<p(— tA) (-Ay" +b(f")) 
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with the entire function (p(z) = (e z — 1 )/z. The matrix-vector product ( p (—xA)v can be approximated 
efficiently by Krylov subspace methods; see lHochbruck & Lubichl(11997l) : ISaadl (119921) . Since the largest 
eigenvalues of the symmetric positive semi-definite matrix A are of size 0(h_ 2 ) in the case of a quasi¬ 
uniform triangulation with mesh size h, the error bounds of Hochbruck & Lubichl 0 997 ) show that 
using m Krylov substeps, each of which requires one multiplication of A with a vector, with a step size 
T = 0(m 2 h 2 /1 loge |) the matrix-vector product <p(— tA)v is approximated with an 0(e) error. For large 
m ~ h~ l the effective step size T e ff = T /m can thus be of size 0(h ), as opposed to only 0(h 2 ) for the 
explicit Euler method or standard explicit Runge-Kutta methods. Runge-Kutta-Chebyshev methods, 
such as those in Abdulle (2002), achieve the same increase in the effective step size. 


5.3 Bulk-surface splitting integrators 

Especially in situations with different time scales in the domain and on the boundary, such as fast 
reaction-diffusion on the surface and slow diffusion in the bulk, it appears attractive to use splitting 
methods that treat the problem in the bulk and on the surface separately. There are basically two ways 
in which such a splitting can be done: 

• Bulk-surface force splitting: Write the stiffness matrix as A = Aq + Ap, where A_q and Aj- are 
the stiffness matrices corresponding to the bilinear forms aq(-,-) and «r(v) defined by the bulk 
and surface integrals, respectively. We use a corresponding decomposition for the load vector 
b(r) = hq (t) +bj-(f). Then a time step with the Strang splitting for this decomposition reads as 
follows: 

1. Over half a time step solve Mu(f) = — Aqu(f) + bq(fo) for 0 ^ t ^ t\n = to + jT with initial 
value u°, to obtain the final value u 1 / 2 ’ - . 

2. Over a full time step solve Mu(f) = —Aqu(f) + bq (t 1 / 2 ) for 0 ^ / < t\ = to + T with initial 
value u 1 / 2 - , to obtain the final value u 1//2,+ . 

3. Over half a time step solve Mu(f) = — Aqu(f) + bj-ffi ) for t f i 2 f t f t\ with initial value 
u 1//2 + , to obtain the final value u 1 . 


With a diagonal mass matrix M, the first and third substep only change nodal values that corre¬ 
spond to surface nodes. In the second substep, bq (t] j 2 ) can be replaced by j (bq (t<f) + bq (t \)). 

• Bulk-surface component splitting: We write the nodal vector as 



where uo and Ui contain the approximate solution values at the inner nodes and boundary nodes, 
respectively. We write similarly 


M = 




f Aoo AoA 
VA10 An ) ’ 



A time step with the Strang splitting for this decomposition reads as follows: 

1. Over half a time step solve Miili(f) = — AnUi(f) — AjoUq +bi(fo) for 0 ^ t < t\j 2 
initial value Uj, to obtain the final value u [ 


with 
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2. Over a full time step solve MoUo(f) = — AooUo(f) — AoiUf/ 2 +b 0 (fi/ 2 ) forO ^ t ^ t\ with 
initial value Uq, to obtain the final value Uq. 

3. Over half a time step solve MiUi(f) = — AnUi(f) — AjoUq + bi(fi) for t\ n ^ f ^ fj with 
initial value u [ /2 , to obtain the hnal value u|. 


In the second substep, bo(?i n) can be replaced by j(bo(fo) +bo(fi)). 

For both types of splittings, the local error (that is, the error after one step starting from the e xact so¬ 
lution) is of order O ( T 3 ) under suffici e nt regularity assumptions on the solution; cf., e.g jDescombes & Schatzman 
( 2002 ): Hansen & Ostermann ( 2009 ): Jahnke & Lubich (12000l) for different techniques to bound the lo¬ 
cal error of the Strang splitting. This yields an 0( T 2 ) global error bound uniformly over bounded time 
intervals in a norm in which the method is stable. 

The force splitting is Li(f2) ©Li(r)-stable: with the symmetric positive definite matrices A_q = 
M-'/^lVr 1 / 2 and A r = M _l / 2 ArM“ 1//2 we have (taking here b(f) = 0 for simplicity), 

M^u 1 = e“5 A r g -|A r m i/ 2 u 0 


and hence 

(u^Mu 1 < (u°) r Mu 0 , 

where we note the /z-uniform equivalence of the norms (see Sec tion l3.1.51 ) 

(u^Mu ) 1 / 2 = \uh\lM ~ |m/i| = {WuhW^a) +Ml|M/ ! ||i 2 ( r)) 1/2 - 

Stability is not evident for the component splitting, since here we have (again for b(f) = 0) a com¬ 
position of exponentials of non-symmetric matrices: 


M'/V =exp(-^f x ° t ° 

1 2 VAio An 


exp —T 




Aoi 

0 


exp 


0 

Aio 


0 

An 


M 1 / 2 u°. 


We have, however, the following stability bound in a stepsize-dependent norm related to the energy 
norm. Here, y 1 = M 1 2 u 1 denotes the mass-scaled numerical solution after one step of the above com¬ 
ponent splitting method with b(f) = 0, and || • ||o denotes the Euclidean norm or its induced matrix 
norm. 


Lemma 5.1 For the matrix 


/(I_ e -rAoo)-l/2 o 0 

V Lio (l-e-^r 1 / 2 ) \ 0 a[{ 2 

with the off-diagonal block Lio = (I — e“ T ^ ll//2 ) 1 / 2 (I + e“ T ^ 11 / 2 )“ 1//2 (A u l '"AioA 00 1 ' /2 ), which has 

IIL10II2 < 1 . we have 

l|Ly 1 ||2<||Ly°|| 2 . 

Proof. For ease of notation we omit the hats on Aoo etc. in the following. We denote the exponential 
matrices 


Eo(t) 

Ei(t) 



^ e -rAoo -(I_ e -^A oo)Aoo l Ao ^ 

(I 0 \ 

V-(I-«" tAlI )AnA,o 
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so that the propagation matrix of the Lie-Trotter splitting is 

Sue =Eo(t)Ei(t) 


and that of the Strang splitting is 

SStrang = Ei (t/2) E 0 (t) E, (t/2) = E, (t/ 2) Sue El (t/2)- 1 . 
It turns out that the block diagonal matrix 

_ Ai-^Aoo) 1/2 Aoo 1/2 0 \ 

0 e TA ‘i/ 2 (I-e- TA ‘i) 1/2 A u 1/2 i 


T = 

V 

symmettizes the Lie mattix: 


S-T-'SueT 

equals the symmetric matrix, with the abbreviation X T = (I — g _,rA oo)V 2 a J^AqiA , 


S = 


e -rAoo + X r (I - e~ tAn )X -x r (l - e -TA u )I/2 e -TA u /2 


-tAi 


l/ 2 (l_ e -TA„ 


)!/ 2 X 


-tAi 


We next show that ||S IU ^ 1- Since the matrix is symmettic, this is equivalent to showing that the nu¬ 
merical range {v T Sv : l|v|| 2 < 1} is contained in the interval [—1,1], With v r = (vf. vf) the calculation 
yields 

v T Sv = Vq e _TAoo v 0 + vJX T (I — e~ xAn )Xv 0 

-2vJX r e“ TAll / 2 (I-e“ TAll ) 1/2 vi+v[e“ TAll Vi, 
where the mixed term is estimated by 

—2vq X r e“ TAll//2 • (I - e- tA " ) 1/2 Vl < vJX r <r TAlI Xv 0 + vf (I - e“ TAl >i, 
so that some terms cancel and we obtain 

v r Sv < Vq e^ TA °°Vo + Vo X r Xvo + vf vq. 

Since the Schur complement Aoo — A 0 | A, f A 2 , is symmetric positive definite, we infer that also the 

transformedmatrix I—Aqo^^oiAjj 1 ^-Aj/^AqjAoo^ 2 is positive definite and hence ||Aqo'^^Aoi Aj j 1 11 2 
1, so that on recalling the definition of X we can bound 

VqX^Xvq < Vo (I —e -TA °°)v 0 . 


This gives us finally 


v r Sv < VqVo + v[v, = II v||| 


and in the same way we also obtain 

v r Sv ^ —||v||l- 

Hence we have shown that ||S II2 ^ 1. Since 


Sstrang=E 1 (T/2)S Li eE 1 (T/2)- 1 =E 1 (T/2)TST- 1 E 1 (T/2)- 1 , 

we obtain with L = (Ej (t/ 2)T) _1 , which is the matrix L given in the statement of the lemma, that for 
y 1 = Sstrangy 0 we have the bound 

l|Ly 1 1| 2 = 1 1 SLy° 11 2 < ||Ly 0 || 2 , 


which concludes the proof. 


□ 
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